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Abstract 

We present an approximation scheme for calculating observables and strength 
functions of finite fermionic systems at finite temperature such as hot nu- 
clei. The approach is formulated within the framework of the Hubbard- 
Stratonovich transformation and goes beyond the static path approximation 
and the RPA by taking into account small amplitude time-dependent fluc- 
tuations around each static value of the auxiliary fields. We show that this 
perturbed static path approach can be used systematically to obtain good 
approximations for observable expectation values and for low moments of the 
strength function. The approximation for the strength function itself, ex- 
tracted by an analytic continuation from the imaginary-time response func- 
tion, is not always reliable, and we discuss the origin of the discrepancies and 
possible improvements. Our results are tested in a solvable many-body model. 
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I. INTRODUCTION 



Mean-field approximations [|T],3, such as the finite temperature Hartree-Fock (HF), are 
standard for describing nuclei at finite temperature. Collective excitations of the system 
are obtained by linearizing the time- dependent HF equations around the static HF solution, 
leading to the finite temperature random phase approximation (RPA) |Q. However, this 
treatment is inadequate for situations where various nuclear configurations compete and 
have comparable free energies, which causes observables to fiuctuate widely about their 
mean values. This is the case, for example, in the vicinity of a shape transition from a 
spherical to a deformed nucleus, where various discontinuities predicted in the mean-field 
approach are smoothed in the finite nuclear system by the presence of fiuctuations. Large 
fiuctuations in the nuclear shape also play an important role in giant dipole resonances 
(GDR) whose frequency is strongly coupled to the quadrupole deformation. In the adiabatic 
approximation, the observed GDR strength function is obtained by integrating the strength 
function that corresponds to each quadrupole configuration over all possible configurations 
weighted by their respective Boltzmann factor [Q. 

It has been shown [0,^ that a systematic description of fiuctuations can be obtained in 
the framework of the auxiliary-field path integral (Hubbard-Stratonovich transformation 0). 
In this framework the original many-body propagator (in imaginary time) is decomposed 
into a superposition of one-body propagators describing non-interacting fermions moving in a 
fiuctuating external time-dependent potential (auxiliary fields). The expectation value of any 
observable is represented as a weighted average of its expectation values in the corresponding 
non-interacting systems. Applying the method of steepest descent one can obtain self- 
consistent mean-field approximations at finite temperature. By taking into account small 
amplitude time- dependent fiuctuations of the auxiliary fields around the mean-field solution 
in the Gaussian approximation one rederives the RPA 

The path integral representation also constitutes a starting point for exact numerical 
solution of the many-body problem where the auxiliary field integration is performed by 
Monte Carlo techniques that were recently developed for strongly correlated electron systems 
P] and for the interacting nuclear shell model M. In the latter case a practical solution 



to the Monte Carlo sign problem [|T0|], which is generic to all fermionic systems, enables 
the study of nuclear properties in medium and heavy mass nuclei using realistic effective 
interactions. The auxiliary field formulation also suggests new approximation schemes which 



are non-perturbative, the simplest of which is the static path approximation (SPA) |TT 



- [0. Here, the path integral is approximated by summing over the time-independent 
fields only, weighted by the appropriate Boltzmann factor. This amounts to averaging 
over all possible static mean-field configurations rather than the self-consistent ones alone. 
The SPA has been used to calculate free energies and level densities in nuclei and was 
found to be superior to the mean-field approximation. In particular, it accounts for the 
enhancement of level density due to thermal fiuctuations of the shape [jl3[ . More recently 



the approximation has also been applied to the calculation of strength functions |T^. At 
high temperatures the SPA partition function approaches the exact result. However, as the 
temperature decreases the SPA becomes inaccurate since time-dependent fiuctuations about 
the mean-field configurations can no longer be neglected. This is manifested especially in 
the strength function where even its first moment is significantly underestimated flGl . 
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Recently, a method to improve the SPA has been proposed for the partition function 
ITSl^ll. In this approach contributions from time- dependent fields in the neighborhood of 
each static field are incorporated perturbatively to the second order in their amplitudes. 
When one considers such time-dependent fluctuations around the equilibrium configuration 
only, one obtains the RPA corrections to the partition function. However, in the proposed 
approach such small amplitude time-dependent fluctuations are taken around each static 
configuration of the auxiliary fields and the static integration is still fully retained. At even 
lower temperatures some of these time-dependent fluctuations can become unstable and 
the approximation breaks down. However, this happens only at temperatures below the 
phase transition, when the temperature drops below the largest imaginary RPA frequency. 
We remark that as a saddle point develops in the nuclear free energy surface below the 
transition temperature, it is still possible that the imaginary RPA frequencies are small 
enough in magnitude that the time-dependent fluctuations are all stable (although there is 
an unstable direction in the static free energy surface). This approximation scheme has been 
applied to free energy and level density calculations in simple models and shown to work 
well down to low temperatures. 

It is interesting to investigate whether incorporating small time-dependent fluctuations 
provides a significant improvement over the SPA evaluation of quantities other than the par- 
tition function. This paper will explore the validity and applicability of this scheme, which 
we term the perturbed static-path approximation (PSPA), for the calculation of expectation 
values of observables at finite temperature as well as of strength functions. Previous work 
employed a formulation of the PSPA based on ordinary quantum-mechanical perturbation 



theory |l9| , p0| , which was specialized for free energy calculations and is not easily extended 
to other quantities. We therefore reformulate it in a general way, using many-body methods, 
and apply it to the calculation of observables and strength functions. By testing the PSPA 
in a solvable model we find that the PSPA results agree closely with the exact solution for 
observable expectation values and for low moments of the strength function. The strength 
function itself is approximated well at high temperatures and also at low temperatures for 
a certain regime of the model's parameters, although the low-temperature results generally 
do not improve much on the SPA. We discuss possible improvements. 

This paper is organized as follows. In Section II we review the auxiliary-field path integral 
formalism and in Section HI we formulate the PSPA for the partition function. Section IV 
discusses the application of the PSPA to the calculation of expectation values of one- and 
two-body operators. Finally, in Section V we present a treatment of the strength function 
and its lowest two moments in this framework. We illustrate the PSPA and compare it to 
other approximations in a simple model whose exact solution is given in Appendix B. 



II. AUXILIARY FIELD PATH INTEGRAL 

In this Section we briefly review the auxiliary-field path integral representation of the 
imaginary-time evolution operator for an Hamiltonian with two-body interactions (Hubbard- 
Stratonovich transformation) [@,^. We consider a system of interacting fermions, and assume 
for simplicity an Hamiltonian of the form 
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H = K--xV' , 
where K and V are one-body operators 



(1) 



(2) 



creation and annihilation operators for a set of single particle states i, and x is a 
coupling constant. Our results can be easily generalized to a superposition of such separable 
interactions 



We remark that an arbitrary two-body interaction 

ijkl ii \ k / ijkl 



(3) 



aja\ai 



X! - H ^cpPaPp , Pa = alttj , pp = alai 



(4) 



«=i 



a/3=l 



can be brought into the form (j^) by diagonalizing Uajs, with the number q of separable 
interactions equals at most to the square of the number of single particle states. In order 
to obtain a path integral representation of U = exp{—f3H) we divide the imaginary-time 
interval [0,/5) into sub-intervals of length e = f3/N and write 



N 



U=l[ e-^^+^^^^ 



n=l 



N 



Y[ e-^^e^'^^^ X (l + 0{e 



.n=l 



(5) 



At each time slice n we use the identity 



(6) 



valid for any operator d with a bounded spectrum and A > (we shall comment on the 
case A < below), to replace ea^^ ^ by an integral over an auxiliary variable C,n- This is the 
Hubbard-Stratonovich (HS) transformation 0, resulting in 



(7) 



In the limit — > oo, becomes a field ^(t = ne) = C,n and we obtain the continuous version 
of the HS transformation 



(8) 
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where is the (imaginary time) propagator for a time-dependent one-body Hamiltonian 
H^ = K-xaT)V 



f/^ = Texp - j dr [K-xa-r) V] 



(9) 



and T denote time ordering. Eq.(^ provides a representation of the many-body evolution 
operator as an average over one-body evolutions which correspond to non-interacting 
particles moving in a fluctuating time- dependent field ^(t), weighted by a Gaussian factor. 
Following it is advantageous to describe the field ^ (r) in terms of its Fourier components, 



^(r = ne) 



(Af-l)/2 

E 

r=-(Ar-l)/2 



(10) 



where (T_j. = a* to keep ^(r) real and uir = ^-nr / [3 are the Matsubara frequencies. Here we 
assume the number of time slices N to be odd. Rewriting the functional integral over the 
field in terms of cxr, we have 



(11) 



where 



U„ = T exp 



V r^o y 



(12) 



The measure in (|lT]) is given by 



N_ dan T-r dalda'l 
' -7= IT 



(13) 



where ar = <7[. +ia'^. The one-body Hamiltonian H^, whose corresponding propagator is Ua 
separates into static and time-dependent parts 



Ha = ho + hi , 

ho = K- xcToV , /ii = -X E ^re'^'-^V . 



(14) 



It is useful to introduce the interaction picture representation U„ of the one-body propagator 
f/o- with respect to the static part of the Hamiltonian /iq 



Ua = e-^'^oTexp 



J drhi ij) 



(15) 



where /ii(r) = e^^^hie We then rewrite (|TT|) as 
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J c/(Toe-^^'^"oe-/3'^o X J V'[a] exp (-^/^E I I') (16) 



V'la] = n ^ ■ (17) 

r>0 ^ 

Eq. ( |T6|) describes the evolution operator for the two-body H as a Gaussian- weighted average 
of one-body evolutions e~^^° corresponding to a static field ctq, multiplied by a correction 
factor which represents the contribution from time-dependent fluctuations of the ^-field 
about (To- In fact, as we demonstrate below, any quantity of interest can be written as an 
integral over its static-field value times a correction factor. The objective of this paper is to 
approximate this factor for various quantities by evaluating the small-amplitude fluctuation 
contribution to the integral over ar{r > 0) and to explore the validity of this approximation. 

For a superposition of separable interactions (|^) one introduces auxiliary-field variables 
cr° corresponding to each but the preceding development remains unchanged. We point 
out that for X < in (|l|) one should use 

e-^'" = ^ J d^e~'^\^''^^' (18) 

— CO 

(A > 0) instead of (1). 

III. PARTITION FUNCTION 

We now consider the partition function Z = Tre~^^. We work in the grand canonical 
ensemble and set the chemical potential = to keep the notation simple. It is convenient 
to choose a ao-dependent basis for the Fock space in which ho is diagonal: 

^0 = J2'^i(^o)aMoyai{ao) , V = J2'>Jiji^o)aii'^oVaj{c^o) ■ (19) 



We can write Z using (0) as 



Z = TtU = J V[a] exp (^-\xPj2 I Trf/. = J V[a 



(20) 



or in a form more convenient for our purpose (using (|TBp): 

Z = Mf daoe-^-^'^koC'o - M I rfaoe-^^o('^-o) . (21) 
\ Ztt J \l Zn J 

Here 

Co = Tre"^'*" = n [l + e-'^^'('"o)] (22) 
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is the partition function corresponding to the static part ho and 



a = J V'[a] exp (^-xf3 ^ | a, i-Tr (e-^'^^W.) = J V'[a] exp l^-xP I (W.)o (23) 

is the correction factor to Co due to the time-dependent fluctuations of the field about (Tq. 
We use the notation (0)o = Tr(e~'^'*°0)/Co to denote the thermal average of an observable 
O with respect to Hq. The effective static-field free energy Fq is defined by (^) to be 

oo) = \x<yl - ^ log Co - ^ log Co ■ (24) 

The representation (^TJ) of Z is a starting point for various approximations. The mean- 
field approximation (MFA) is obtained when the contribution of the time-dependent paths 
to Z is neglected by setting hi = 0, implying Co = I5 and the integration over the static 
fields (To is performed in the method of steepest descent. This amounts to approximating 
the path integral by the contributions of the static paths ctq that minimize the free energy 
Fq in (|2^). It is easy to show 0,^ that this minimization condition is 

^0 = {V)o = J2^iii^o)f^i^o) , (25) 



where fi = (l + e^"*)"^ are the Fermi occupation numbers, and that the solution ctq of ( p5D is 
the Hartree mean field. One can improve on the MFA result by performing the integration 
over ar in the expression (|23|) for Co also by steepest descent. The saddle point is now 
given by o"o, (Jr = and one obtains the finite temperature RPA corrections to the partition 
function 0. 

The static-path approximation (SPA) [Tl| - [|l7l is obtained by again setting hi = 
but now the integration over o"o in (|2T| ) is performed exactly, thus approximating the path 
integral by the contributions from all static paths: 




rfaoe-^^^'^oCo = \ — I rffxoe-^^^"" H + e-''^'^"")] . (26) 



The SPA is expected to become exact at high temperatures since one can use the one 
time slice approximation in (|^ with an error of 0{j3'^) that vanishes as T — 00. This 
method is advantageous to the MFA since it takes into account exactly large amplitude static 
fluctuations around the mean field. However, it neglects the time-dependent fluctuations 
which constitute the RPA corrections. This shortcoming of the SPA can be remedied if the 
exact integration over the static paths is supplemented by evaluation of Co in the limit of 
small amplitude time-dependent fluctuations about each static value (Tq. This scheme, the 
perturbed static-path approximation (PSPA), is the focus of our work. It was introduced in 
Refs. [0 - 1^ for the partition function Z. Our approach is different and has the advantage 



that it can be generalized to the calculation of observables and response functions as shown 
in the following sections. We first illustrate our method by rederiving the corresponding 
expression for the partition function. We expand log(Wo-)o in ( ^31) 
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1 ^ 

+ E <'r<y sj drdr'e'-^-e^-^-' [(TV{t)V{t')), - (K(r)),(V^(r')).] 



+ O {a') , (27) 

with V{t) = e'^^''Ve~'^'^° in the interaction picture. The time-ordered averages are calculated 
using the finite-temperature Wick's theorem where the subscript c means that only 
the connected diagrams should be summed up in the diagrammatic representation of this 
expansion. 

The first term in (^) vanishes since {V{t))c is r-independent. For the second term 
Wick's theorem gives 

{TV{t)V{t'))^ = - E v,,v,,g^{r' - r)g^{r - r') , (28) 



where the unperturbed temperature Green's function is given by 

1 oo p-iUkir-r') 

gKr - r') = -(a.(r)al(r'))o = E • (29) 

fc=-oo '^'^'^ 



with frequencies z/fc = (27r + l)k/ [3. Using (^8| - p9|) we obtain for the double integral in (p7|) 



oo 2 1 
drdr' ■■■ = -5r-s E E ~ 7 ■ — T ■ (^0) 

^ ^r^^ ^Vk - ei IVk - (Cj - lUJr) 



The infinite sum over k in ( pOf ) is calculated using the frequency summation technique . 
The essence of this method lies in the observation that the points z = ivj. are the poles 
of the function —f3/{e^^ + 1) with residue of one and the sum therefore equals the contour 
integral 



oo 



E 7 7 — / n ^ 7 T ; (^1) 

where C encircles the imaginary axis. This contour can be continuously transformed into a 
circle centered at the origin of an arbitrarily large radius, which is deformed at two places 
to include the poles at 2; = and z = ej — iur- The residue theorem then gives 

/ drdr' ■ ■ ■ = —dSr E VuVa—-^ — , (32) 



with Ajj = ej — ej. Using (|32D in (^) we have 



(W^)o = exp -Xl3 E «r I (Tr ? 

\ r>0 / 



(33) 
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with 



«r((Xo) = X E ^ij^n /' J-,, = X E ^ipji ^\2 ^l]^'^ ■ (34) 
The correction factor Q in (p3D is now given by a Gaussian integral and we obtain 

c; . n (1 . ^ n ^j:^ - ^y-i^ . (35) 



u 



The second equahty in (^) defines the frequencies r2!y(/?,(Jo) through |20 



where the prime in 0' cind E' restricts the product or sum to pairs {i,j) that satisfy i < j 

ij ij 

and Aij ^ 0. Note that there are as many Qi, in the numerator of (^) as there are Ajj 
in the denominator. The third equahty in (pSj) uses the infinite product representation 
sinhx = X n (1 + x^/vr^r^). Together with (|2T|) and ( ^2]) we finally have 



r>0 



UnJ -^M ^ n^sinh^ ^ ^ 

This is a closed-form expression which corresponds to the limit N ^ 00 {N is the number of 
imaginary-time slices, see (|^)) since the infinite product in (|35|) has been performed exactly. 
Hence the PSPA result does not contain errors originating from a discretization of [0,(3) 
into sub-intervals of a finite length, as is the case in a Monte Carlo evaluation of the path 
integral. 



LU = ±^ly are the roots of 1 -|- E VijVji{fi — fj)/ {^ij + iuj) = 0. It can be shown |^ that 

ij 

for the value ctq = (Jo{P) which minimizes Fo(/5;cro) (see Eqs. (^) and (p5|)), these roots 
ifij, are the RPA frequencies, i.e. the frequencies of small amplitude oscillations around 
the equilibrium configuration. For an arbitrary static ctq, the frequencies u = iiQu solve the 
generalized finite temperature RPA equations obtained by replacing the equilibrium config- 
uration by the arbitrary ao (see Appendix A) . We note that the Gaussian approximation in 
(p3|) leads to a convergent integral only if 1 -|- > for all r, i.e. < 00^ for all r and u. 
When all RPA frequencies are real, these conditions are always met. In particular this is the 
case if ctq is a local minimum of the static free energy i.e. a stable mean-field configuration 
(Jo (Thouless theorem |^). At zero temperature Ur = and an imaginary RPA frequency 
would lead to a breakdown of the Gaussian approximation. However, at finite temperature 
imaginary RPA frequencies do not necessarily lead to instability. If the largest modulus 
of all imaginary RPA frequencies is below Ui = 2ttT, then the quadratic fiuctuations in 
(Tr are still stable. An instability in the ai direction occurs when the magnitude of one of 
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the imaginary RPA frequencies crosses 2'kT. The RPA frequencies depend both on (3 and 
(Jo. In particular, for any temperature there exists a static field do = o"o(/5) such that 



-CUT 



-{2tt/ f3Y for some causing the correction factor C^'q to diverge. Note 



that at T = instability occurs as soon as some Vty becomes imaginary whereas at T > an 
instability develops only for a sufficiently imaginary Vt^. These instabilities can be ignored 
when they occur at fields for which the static free energy is large. In practical applications 
the PSPA breaks down only at temperatures that are significantly below the shape transi- 
tion temperature, when the saddle point in the static free energy surface becomes unstable 
to small amplitude time-dependent fiuctuations. 

In the case (^) where the interaction is a sum of several separable interactions we have 
an auxiliary- field variable cr" for each Va. Assuming Xa = X (this is always possible by 
redefining Va), expression (p3|) becomes 



exp -X/? ^ ^r' 

\ r Q7 



(38) 



with 



07 



>o) =xYl 



1-j 



V V ■■ — 

11 J* A I • 



The correction factor Q'q is then generalized from (|35D into 

C= ndet(l + a,)-^ 

r>0 



(39) 



(40) 



where is the g-dimensional matrix defined in (pOf). As shown in Appendix A 

det(l + ttr 



(41) 



where VLy are again the RPA frequencies at finite temperature. We then obtain for the PSPA 
partition function an expression similar to (|37|), except that the integral over ctq is replaced 
by an integration over q static fields cTq . 

When one or more of the Xa in (0) are negative and a representation of the type ([T8| ) is 
used in the HS transformation, the one-body Hamiltonian in ([l9|) becomes non-hermitean, 
so that its eigenvalues ei(cro) are in general complex and the associated one-body partition 
function can be negative. It was recently pointed out [0] that the static fields which 
correspond to these "repulsive" terms in the interaction do not represent large amplitude 
thermal fiuctuations, and can simply be treated in a saddle point approximation. This 
amounts to keeping the exact integration over the static "attractive" fields but taking the 
mean-field solution for the static "repulsive" fields for every configuration of the "attractive" 
fields in the integrand. 

An alternative way to approach "repulsive" interactions is to subtract from the Hamilto- 
nian a term proportional to A^^ (A^ is the particle number operator) with a coefficient large 
enough so that the two-body part of the Hamiltonian becomes a negative-definite quadratic 
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form (as in the case where all Xa in (U) are positive). However, by doing this we also increase 
the magnitude of "attractive" terms and it remains to be investigated whether the SPA and 
PSPA would still work well for such modified Hamiltonians. Note that the subtracted term 
is just a constant for a fixed number of particles. 

To illustrate our results we apply the formalism to a simple many-body model (a variant 



of the Lipkin model |2^) based on a U(2) algebra which is defined and solved in Appendix 



B. At each cxo we have 

/lo = 2eJ, - 2xaoJ. , h = -2% ^ a.e^'^'-V, . (42) 

Thus the single-particle Hamiltonian corresponding to Hq has two (7-fold degenerate levels 



ei = ±e where e = ye^ + x^'^o; so Aij = 0, ±2e. The matrix corresponding to = 2Jr^ is 
block-diagonal with g 2 x 2-blocks 



Vi 



( -sm24> cos 20 \ . X(To e , . 

jn — \ „ T • o r 5 sm 20 = , cos 20 = - (43) 

^ y cos20 sm20 y ' ^ g ' e ^ ^ 

and the matrix 

Pi=( "^^"^ ''""^1 (44) 
I — sin cos y 

diagonalizes the 2 x 2-blocks of ho. For ^0 and Co we obtain 

C.= (2cosh-) . C = ^-i^. (45) 

where = 4e^ — (4(7e^x/^) tanh(/?e/2), so there are only two RPA frequencies ±Q. The 
saddle-point mean-field equation (|25|) is 

tanh ^ = — (46) 

and has a single spherical stable solution (Tq = for /? < (3^, which becomes unstable for 
P > Pc and bifurcates into two deformed minima ao = ±cto(/3). Here Pc is the transition 
temperature 

/5c = - log , K = — . (47) 

e K — 1 e 

The dimensionless parameter k thus determines the mean-field behavior of the model. For 
K < 1 there is no phase transition and Uo = is the only solution at all temperatures. 

In Fig. |T] we present the MFA, SPA and PSPA results for the free energy F{P) = 
—P~^ log Z in our U{2) model using e = 1 and g = 10 and compare them against the exact 
result. We consider three cases characterized by different values of k. For ft = 0.5 there is 
no mean-field transition whereas for n = 1.5 and k = 3.0 the signature of a transition at 
temperatures Pc = 1-61 and Pc = 0.693, respectively, is seen in the mean-field curve. The 
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SPA overestimates the exact result especially at low T whereas the PSPA is quite accurate 
at all temperatures. All three approximations converge to the exact result as T — oo. We 
find that for /5 sufficiently small the divergence of Co occurs only at a very large ctq = o"o(/5), 
where the other factors in the integrand of (|37D are vanishingly small, and therefore does not 
constitute a practical problem. The breakdown of the PSPA occurs when (3 becomes large 
enough such that cro(/3) ~ 0. This happens only at temperatures far below the transition 
temperature Tc = 1/ (3c. 

In order to study the effect of the time-dependent ffuctuations we compare in Fig. the 
effective free energy Fo(/?; ctq) in ( P^ ) as a function of the static field (Tq with the SPA free 
energy {C'q = 1 in (p^)) at different temperatures. We consider the case n = 1.5. Above 
the mean-field transition temperature f3c = 1.61, Co ~ 1 the two approximations yield 
similar results. As the temperature is lowered, time-dependent fiuctuations deepen the free- 
energy minimum at ctq = 0. Below these fiuctuations also lower the barrier between 
the two mean-field configurations ctq = ±cto. Thus the PSPA result for the free energy 
improves significantly the SPA result especially below the transition temperature (see the 
middle panel of Fig. |1]). The MFA result for the free energy is the most inaccurate 

especially near the transition where the free energy has comparable values over a broad 
range of configurations a^. 

We remark that in spite of its simplicity, the U{2) model has features that are generic 
to more realistic nuclear interactions (e.g. quadrupole interaction). In particular, the phase 
transition in its mean-field theory is analogous to the shape transitions from deformed to 



spherical shapes that occur in deformed nuclei ||2^. We therefore expect the PSPA to 
improve significantly on the SPA also for more realistic nuclear shell model interactions. 
For interactions that include both pairing and multipole components (e.g. pairing plus 
quadrupole model), the SPA works better when a mixed pairing-density decomposition is 
used rather than just a pure density decomposition 0. A mixed decomposition is thus 
preferable when the PSPA is applied to such interactions. 



IV. THERMAL EXPECTATION VALUES OF OBSERVABLES 

In this Section we consider an observable O and treat its expectation value at finite 
temperature (O) = TT{e~^^ O) /TT{e~^^) in the framework of the PSPA. Although we discuss 
observables of the form O = D and O = for a one-body operator D, any n-body O can 
be treated in a similar fashion. 



A. One-Body Observables 



In constructing an auxiliary-field path integral representation for {D) we can use Eqs. 
P) and (pOl) to get 



Tr(f/,Z}) 



Alternatively, to facilitate our calculations we can use (p!6D and (21) to write 



(48) 
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I daoe 



where 



Do = ^ f V'[a] exp f-x/^E I A {U.D)o . 



(49) 



(50) 



The effective static-field free energy Fo(/3;(To), defined in (0), has been calculated in the 
previous Section (see Eqs. (^) and (|35|) ). Our aim here is to calculate Dq. In a diagrammatic 
representation of the perturbation series for (lAcjD)^ the connected diagrams can be easily 
shown to factor out at each order 



Thus we can write Dq itself as a weighted average of {UaD)^ 



(51) 



Da 



JV'[a] exp i-xPEl CTr ? {UM^^D), 

\ r>0 / 

/I)'[(T]expf-x/3E \<yrA (W.)o 

V r>0 / 



(52) 



(Wo-)o is given in (^) and (0), whereas for (U^D)^ we expand 

1 r 

{U,D), = {D)a + -X^ E ^r(Tsj dTdT'e'^^^e''^^^'{i:V{T)V{T')Dm, + 0{a^) . (53) 

Notice the absence of first-order terms in cr^ since they vanish upon the integration in (p^. 

Before continuing we introduce the following notation. Using ( p3D and (^) we define a 
matrix A(cro) by 



Co = / ^^'Mexp 



r>0 



(54) 



where the vector a has — 1 components (crj, (72, ■ ■ • , O"" , cr'-^, ■ ■ ■) and A is diagonal with 
elements (1 + ai, 1 + 02, ■ ■ ■ , 1 + ai, 1 + 02, ■ ■ ■) with given in (|3^). We also use (|53D to 
define &(cro) and a matrix -B(cro) by 



Using this notation we find 



(SPA) 



(55) 



(P5PA) 




(56) 
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Recall that at each do we choose a different basis ([I9D for the Fock space in which we 
write 



The zeroth-order term is then 

6(0-0) = (^)o = ^diifi . 



(57) 



(58) 



In order to obtain an expression for Tt{A~^B) in (^) we calculate the double integral in 
(^) using Wick's theorem and the frequency summation technique |T]J^. Notice that since 
A is diagonal it is sufficient to consider the case s = —r. The calculation is similar to the 
one performed in the previous Section and yields the result 



ijk Ty^O 



(59) 



where 



jijk 



fk 



In the cases ej = etc. it is understood that I^^^ = lim I^^^. 

6 q *■ e i 



(60) 



The sum over r 7^ in ( ^Qf ) which consists of A^ — 1 elements can be expressed in a closed 
form for ^ 00 by applying once more the frequency summation technique. Writing 
= a~^{iujr) and J*-'^ = P^^{iUr), we observe that the sum is carried out over the points 
z = iujr which are the poles of the function (3/{e^^ — 1) with residue of one. Therefore 



(61) 



where the contour C encircles the imaginary axis. To evaluate the integral we transform C 
into an arbitrarily large circle, deformed to include the poles ai z = Aij, z = —Ajk, and 
z = ±i7,, for all and obtain 



(SPA) 




Here 



(PSPA) 




^ ] diifi ^ y '^ij'^jkdki 



2(3 



ijk 



E 



/3 



Pijk 



PijkjO) 

1 U (Za- Za') ' n(-^a') 
a'j^a a' 



+ 



(62) 



Pijk{Z) 



^ik 

Aiji Aj/C5 i^^u 



A 

A. 



ik 



[z + A,-fc) + f, 



n'(-'-AL) 



Im 



(63) 
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and the prime in ' restricts the product to pairs (/, m) satisfying / < m and Aim 7^ 0. 

Im 

In the case of g > 1 separable interactions in @ the calculation is more complicated 
since the matrix A{ao) in ( |5BD is block-diagonal with blocks of dimension 2q (see Eqs. (|5BD- 
(^)). However, the general form (^) including the summation over r still holds and the 
subsequent use of frequency summations to get a closed-form result can be carried out in 
this case as well. 

We test the approximations (|62|) by applying them to the calculation of (J^) in our 
U{2) model {{Jx) = (Jy) = both exactly and in the above approximations). The matrix 
corresponding to D = is block-diagonal with g 2 x 2-blocks 

d., = -( sin 20 \ 

2 I sm 2(j) - cos 20 / ^ ^ 



in the notation of (13)- (44). Thus 



(J.)(^^-^) = ^/ff / daoC'^^^^o (2cosh^y' X ^cos20tanh^ . (65) 

The more complicated expression for {Jz)^^^^^^ is similarly obtained using (^)-(|45|) and 
(|6^)- (|6^) . The SPA and PSPA results are presented in Fig. ^ together with the mean- 
field calculation derived from a steepest-descent evaluation of (§5|). Above the transition 
temperature all approximations agree fairly well with the exact result (p8|). For lower 
temperatures the MFA and SPA significantly overestimate the exact result whereas the 
PSPA works well down to the breakdown temperature. 



B. Two-Body Observables 

For simplicity we assume a two-body observable of the form O = D'^D where D is a 
one-body operator. The calculation of the expectation value of such an observable, albeit 
more complicated, can be carried out along the same lines. As for {D) we have 

where 

0, = ^ I V'[a] exp {-xP E I <'r ?] {WD^D)o . (67) 

Following similar steps we get an expression of the form (|56|) with the same A{(Tq) but 
different B{aQ),h{aQ). After further manipulations we obtain the final result 

Or^' = h{a,) , 

Or''^) = %o) + ^Et.(ao), (68) 

u=l 
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where 



^(^o) = J2 [4jdjifii^ - fj) + 4idjjfifj 

ij 



ijkl 



ijkl 



ijkl 



t4(cro) = -^VijVjkdkid\jiSfj 



(1) 

ijk 1 



ijkl 



ijkl 

- '"ij'"kld]kdHS-jli . 
ijkl 



^6(cro) 



(69) 



sijL in (H) are given by 



<lljk^) = Vijk{z) of (63) , 

qiU^ = -{h - mfk - fim\z' - Aj^) , 



V = 1,2 



pq 



Ajj? ^jki i^j/ • 



(70) 



To test the approximations ( |68D we apply them to the calculation of (Jj), {Jy) and (J^) 
in our U{2) model. The matrices corresponding to D = J^, Jy are block-diagonal with g 
2 X 2-blocks 



(j^) _ 1 / —sin 20 cos 20 
2 I cos 20 sin 20 



(J,) _ 1 / - z 



(71) 



whereas d^/j''^ is given in (^4]). We therefore have 



2\{SPA) 



1 /x^ 



I tiaoe-^^^'^o cosh x {D^^ 



2\(SPA) 



(^^) 



2n(5pa) _ g 



f3e , ^7 



/5e 



sin 20 tanh 1 — cos 40 tanh hi 



4 " 2 



= ! tanh^ ^ + 1 



8 V 



2\(SPA) _ 9 



Pe , g 



(31 



4 



cos 20 tanh 1 — —cos 40 tanh hi 



(72) 
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Expressions for I^D'^)^^^^^^ can similarly be derived using (^3])-(^) and (|68| - [7TD . These 
results are presented in Fig. ^ together with mean-field calculations for different values of 
K. The MFA and SPA exhibit large deviations from the exact results in most cases even 
above the transition temperature and when no transition occurs. The PSPA shows the 
best agreement and in the only case where its deviation is appreciable {{Jy) at k = 1.5), the 
other approximations give qualitatively wrong results. Note that the quantitative differences 
between the various approximations depend on the observable being calculated. In the no- 
transition case K = 0.5 all approximations agree for (J^) (as they do for (Jz)) but not for 

The expectation value of the Hamiltonian H = 2e — 2xJ^ ( piP itself is shown in Fig. 
H (note that the additivity of the average, (Oi + O2) = {Oi) + (O2), is preserved in MFA, 
SPA and PSPA). The PSPA result is still superior but the agreement of the SPA with the 
exact result is also quite good, contrary to what it predicts for (J^) and (J^) when taken 
separately. 



V. STRENGTH FUNCTION 

A. Linear Response Theory 

In this Section we are interested in the response of the system to an external perturbation 
D as reflected by its thermal strength function 

Gi^) = ^ E I I ^ I I' e-^^™5(a; - E„ + , (73) 

mn 

where | n) and £"„ are the many-body eigenstates and corresponding energies. G{uj) is the 
Fourier transform of the real-time response function G{t — t') which originates in linear- 
response theory |l|J^ 

dt 



G{u) = f ^e-*G(t) 
J Ztx 



G{t - t') = -z{TDUt)DH{t')) = --Tr [e-^''TDUt)DH{t')] , (74) 

where the subscript H refers to the Heisenberg picture Dnit) = e*^*De~*^*. If an equi- 
librated system with a Hamiltonian H is perturbed by some external potential V{t) at 
t = to, the change in the thermal expectation value of an observable D with respect to the 
equilibrium situation is 

. t 

6{D){t) = -'-jdt'{[DH{t),VH{t')])Q{t-U) . (75) 

to 

Considering an external perturbation of the form V{t) = a{t)D = Y^a{t)dijalaj we can 

ij 

write ([75| ) as 
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1 °° 

S{D{t)) = -j dt'G'^it - t')a{t') , (76) 

to 

where 

G\t-t') = -z{[DHit),DH{t')])e{t~t') (77) 

is the retarded real-time response function. The Fourier transform of the latter is related to 
that of (5(D) (t) by 

oo 

6{D){u) = '^a{uj)G'^{u) , G^(u;) = ^ / dte'"*G^(t) . (78) 

a ZTT J 

— oo 

The strength function can also be obtained directly from G^{uj): 

G{u) = -l^-^lmG^{u) . (79) 

This suggests that in order to approximate G{uj) we should approximate G{t) or G^{t) and 
Fourier-transform the result. This has been done in |jl6| where G{uj) has been obtained from 
a static-path calculation of G{t). However, it is very difficult to incorporate the contribu- 
tion of time-dependent paths into the real-time formalism. In contrast, the imaginary-time 
framework discussed in the previous Sections is quite suitable for this task. Thus instead of 
G(t) we consider the imaginary-time response function 

g{T - r') = -{TDUr)DH{T')) = -^Tr [e-^^TDUr)DH{T')] , (80) 

where Dh{t) = e'^^De^'^^ . Q is periodic with period /3 {Q{t + jS) = Q{t)), hence we are 
interested in its Fourier coefficients 

dre'^^-^gir) , (81) 



where Un = 27m/ f3 are the Matsubara frequencies G^{uo) and Qn are related by an ana- 
lytic continuation since their Lehmann representations are determined by the same weight 
function p{uj). Defining 

Viz) = I dcoP^ , 



z — u 



p{uj) = 7 (l - e-^") Y.\(^\D\m)\^ e-^^-5(^ - E„ + E^) , (82) 

mn 

it is easily verified that 

G'^iu) = r(cU + Z7]) , = T{zUn) (83) 
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as —> O"*". We can therefore calculate G^{uj) by obtaining an expression for Qn with an 
explicit dependence on iun, then perform the analytic continuation by formally replacing 
iun ^ u + irj. We can then calculate the strength function G^u) from (|79|). In general, an 
analytic continuation of Qn is not unique since it is based on extrapolating from a discrete 

set of points iUn to a continuum u. The sum rule / duopiuo) = (D'^D), however, selects a 

— oo 

continuation that satisfies T{z) ~ {D'^D)/z as z —>■ oo which is unique In this Section 

we present an auxiliary-field path integral treatment of Q{iu!n) = Qn from which we extract 
several approximations for G{ijj). 



B. Strength Function in the PSPA 

In the auxiliary-field path integral representation Eqs. (|80|) and ([81| ) can be written in 
the form 



Q{iUn) 



Jdaoe-f^^o 



where 



^0-' \ r>0 ) i 



• 



(84) 



(85) 



The time-dependence in D is understood to be in the interaction picture with respect to /iq, 
D{t) = e^^° De~'^^° . Our aim is to calculate QQ{iujn)- The connected diagrams factor out as 
in fl]) ii: 

dTe'^-^TU^D\T)DiO)), = {U^)o J dre'^-^ [{i:U„D\r)D{Q,)),+ {1:U^D\t)),{U,D{Q)), 



= {U^)o J dTe'--^D\T)Dm, + -{U„)^ E ^rCTsj dTdT'dT"e''^-^e'^^-^e'^^^ 
'"^^o 

X ;(Ty(r')\/(r")/^^(r)Z^(0)), + 2(Ty(r')/^^(r)),(TV(r")/^(0))e] , (86) 
where terms proportional to 5„,o, such as 



rfre*'^""(Dt(r))o(/^(0)) 



; 



(87) 



are omitted since they do not contribute to the analytic continuation. Similarly to (^) we 
can rewrite Qo{iuJn) in the form 



Qo{iu)n) = T7 f 1^'W]e 



1 



K^^n) + ^Tr 



(88) 
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where A{ao) is defined in ([5^ ) and b{ao,iu!n), B{ao,iLJn) are given implicitly by (|86D. The 
calculation of these quantities using Wick's theorem and the frequency summation technique 
is similar to those discussed in previous Sections but much more cumbersome and we give 
here the final results. We have 

PSPA X ^ 

u=l 



where 



_ n' 

J i J j J k J I 



. I • \ o jt J fj PI 

ti (do, «u;„) = -2 2^ VijdljiVkid, 

ijkl 

t2{ao,iuJn) = -2^VjkVkidiidljSl]li , 

ijkl 

t3{ao,iuJn) = -'^^VijVjkdlidiiS^^li , 

ijkl 

U{ao,iuJn) = -2^VijVkidljdiiSljli , (90) 

ijkl 



with 



0(2) _ „-l rijkl 

'-'ijkl ~ "r -'n-rr 

0(3) ^ \^ ^-1 jijkl 

'~-^ijkl r ^—rrn 



ijkl / , ^—rnr ) 
ijkl fi 

( Ajj + iut) {Aik + iujt+s) {Ail + i(^t+r+s) 



fj _|_ fk 



{Aij + iujt) {Ajk + iuJs) {Aji + iujr+s) {Aik + iuJt+s) {Ajk + iuJs){Aki + 



-\ ll . (91) 

{Ail + iuJt+r+s){Ajl + iuJr+s){Akl + lUJr) 

The apparent divergences of Iljli, e.g. when both Aji = and r + s = 0, are handled in the 
usual way by taking the limit ej e/. 

We use frequency summations to bring the infinite sums in S^^jli into closed forms: 

5,-.(ao, - - L e^^. - 1 R (.. - ..0 " Y^) ' ^''^ 

a'^a a' 
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where 



pqJ 



pq 



X 



fi 



{Z - Z2){Z - Zs) H —{Z - Zx){z - 2:3) + /fc {z - zx){z - Z2) 



Za = ^ik + iuJn, ^jk, -^ki, , a = 1, 2, 



(93) 



pgL(^)=nv-A: 



pql 



pq 



X 



'Z - Z2){Z - Zg) + -^{Z - Zi){z - Zg) + fj —{z - Zi){z - Z2) 



V1V2 



Za = ^ij, -^jk, -^jl - i^n, -^^u , « = 1, 2, ■ ■ • 

rii = Ajfc, r]2 = \i + ii^n, Vs = ^ki + i^^n , 



(94) 



pgL(^) = nv-A;j 



pq 



X 



— {Z - Z3){z - Z4) + —{z - Z2){Z - Zi) - —{Z - Zi){z - Z^) - —{z - Zi){z - Z2) 
Vl Tl2 Tl2 Vl 

Za = ^iji \k + "^^m ~ '^^ni ~^kh ; « = 1; 2, ■ ■ ■ , 

T]i = Ail + i^n, m = Ajk+ iujn . (95) 

Eqs. (p9|)-(p0|) and (p2D-(p5D constitute our final result for ^o(^'^n)- The analytic continuation 
(^) followed by employing the relation (^9]) result in 



TT 1 



-/3a; 



uo + irf) . 



(96) 



We remark that the use of frequency summations to convert the infinite sums in S'j-"^/ ( |9lD 



into the finite expressions (^) is essential to the extraction of the strength function, since 
an analytic continuation of the truncated sums would result in a wrong functional form of 

The SPA expression has the simple form 



(SPA) 
' 



(97) 



This result illustrates a significant limitation of this approximation, namely that only tran- 
sitions at frequencies u corresponding to single-particle energy differences Ajj(cro) can be 
described. This shortcoming becomes evident upon performing the integration over the 
static field fB4D which results in 



^(SPA)^ 



1 



EE 

ij A 



d'ijdji{fi fj) 
dAij/dao I 



(98) 



0-0=0-0 
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where the second sum is over the values ctq = (Tq satisfying Aij{aQ) = lo. In the case of our 

U{2) model Aij = 0, ±2^Je'^ + x^^l-, hence for | |< 2e the strength function G'~^^'^\ijj) = 
and the transitions at this tu-range are not reflected. 

It is interesting to compare the shell-model Monte Carlo (SMMC) methods with the 
PSPA. While in the SMMC the auxiliary-field path integral is evaluated exactly (except 
for statistical errors), the problem of extracting the strength function from the imaginary- 
time response function is quite difficult due to statistical noise. The strength function is 
calculated in the Monte Carlo using a maximal entropy reconstruction method , but this 
method work well only in some cases. In contrast, the PSPA strength function is extracted 
by exact analytic continuation (although only within the approximation). An additional 
advantage of the PSPA is that the infinite-discretization limit of the imaginary-time interval 
[0, (3) is taken exactly, whereas in the SMMC it is necessary to extrapolate the finite time 
step to zero. The validity of the PSPA for strength functions is tested in Sections V.C and 
V.D. 

C. Moments of the Strength Function 

The moments of the strength function 

oo 

M„= J duj uj''G{uj) , n = 0,1,2,... (99) 

— oo 

provide another measure of the quality of the approximations we develop. Rather than 
integrating over the expressions we have for G{uj), it is more convenient to obtain the 
moments directly in terms of H and D. From G{uj) being the Fourier transform of G{t) it 
follows that 

fl"'G 

M„ = z"^U, (100) 
and differentiating the definition of G{t) (second Eq. in (fTiD) we get for the lowest moments 

Mo = (D^D) , 

M, = ^{[d\[H,D]\) . (101) 

Thus the zeroth moment (total strength) is simply a two-body expectation value, which was 
discussed in the previous Section and found to be well reproduced by the PSPA (contrary 
to the SPA) in our U{2) model (see Fig. |). 

The first moment Mi (Mi/Mq is average transition energy) is also given by the expec- 
tation value of a two-body operator. Taking D = in our U{2) model we can exploit the 
angular-momentum commutation relations to get Mi = 2ie{JxJy)- The latter is calculated 
by a generalization of the results for {D'^D) obtained in the previous Section to the case 
{D\D2): whereas dij in (pO|) is the matrix corresponding to D2, d\j should be taken to be 
that of Di. The results for Mi are shown in Fig. for different values of the mean- field 
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parameter n in p7|). The PSPA agrees well with the exact result. The SPA is good at high 
temperatures but worsens appreciably near the mean- field transition and below it, except for 
the no-transition case {k = 0.5) where it remains a good approximation also at low tempera- 
tures. This comparison is interesting in view of the results of Ref. where two versions of 



the static approximation for the strength function were tested. The first version (called the 
adiabatic approximation) is identical to our SPA. The second version (called the static-path 
approximation in Ref. |TB[) consists of estimating the path integral representation of Q{t) 
in (|80|) by step function paths C,{t') which have a discontinuity at r' = r, 

eM=l!°' (102) 



rather than by the constant ones (uo = po) alone. We shall discuss this approximation 
further below. The PSPA includes the effect of discontinuous paths only to the second order 
in (Jo — Po through the correction factor from small oscillations about the average value 



o-qt/P + (po — cro)(l — t/(3). It was found in Ref. |jT6[ that although the inclusion of the 
discontinuous paths ( |102| ) provided an improvement over the SPA, a good agreement with 
the exact results was not achieved, contrary to the situation in the PSPA case. This suggests 
that the important contribution to the moments of the strength function beyond the SPA 
comes from small-amplitude oscillatory paths rather than from step function paths with a 
large discontinuity. 



D. Strength Function in a Simple Model: Results and Discussion 

We test the PSPA for G{uj) in our U{2) model with perturbing operators D = J^. and 
D = Jy. The results are presented in Fig. |^ and Fig. H, respectively. We consider the same 
three cases studied previously characterized by k = 0.5,1.5,3.0 at different temperatures, 
chosen to be at the mean-field transition in each well as above and below it (a 

transition does not occur for k, = 0.5). Note that the exact result (p9|) consists of a sum of 
(5-functions in the limit t] —>■ 0^ and is therefore singular, as is the case for the MFA result. 
This is not the situation in the SPA and PSPA which involve an integration over the static 
field (Tq. In order to facilitate a meaningful comparison we keep rj small but finite both in 
the exact result, which becomes a sum of Lorentzians of width r], and in the approximations. 
Finally, we plot G{uj) only for positive u since G{—uj) = e~^^G{uj) for an Hermitean D. 

The shortcoming of the SPA expressions (p7D-(p8D is manifested clearly in Figs. |^ and p. 
Since we use e = 1 the SPA strength function cannot refiect transitions with uj <2. Hence 
it vanishes in the range uj < 2 even if most of the strength is concentrated there, as is the 
case for n = 0.5, 1.5 at /5 = 1.7,3.0 (where the SPA result has a shifted peak near u = 2 
to the right of the exact peak). Furthermore, since the PSPA result consists of an additive 
correction to the SPA expression (see Eg. (^9])), this shifted SPA peak leaves its trace in the 
PSPA strength function. For k = 0.5, for instance, even though the main PSPA peak is in 
excellent agreement with the exact one, it is accompanies by a small additional (false) peak 
to its right, left over from the SPA, which becomes larger as the temperature decreases. 
In general, the PSPA works quite well for small k but decreases in quality and becomes 
comparable with the SPA as k increases or the temperature decreases. The MFA generates 
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a sharp peak located near the middle of the broader SPA peak, consistently with its origin 
as a steepest-descent approximation of the SPA integral. 

We mentioned above that in Ref. |T6[ G{uj) was calculated in a modified static-path 
approach which included constant paths with a discontinuity at r in the path integral rep- 
resentation of G{t). It is interesting to note that the resulting strength function in [jl6[ is 



accurate for large values of k and deteriorates as n decreases, in contrast with the PSPA 
result. This suggests that for small k the major contribution to the strength function be- 
yond static fields comes from small oscillations about them. However, for large k these 
small oscillations are negligible and large imaginary-time discontinuities in the static fields 
become important. It is therefore desirable to have an approximation scheme which takes 
both contributions into account. In the following we discuss an approach to this problem in 
the imaginary-time framework and the difficulties it encounters. 

We start from the auxiliary-field path integral representation of the imaginary-time re- 
sponse function 



1 
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X Tr i 



Texp 



- j dr' {K - ^{t')V) 



exp 



2% J dreir) 





T 

J dr' [K - ar')V) 



D 



(103) 



However, rather than use the Fourier decomposition ([T0| ) of ^{t') over the entire interval 
[0,(3) to obtain the form (P^, from which the SPA and PSPA are derived, we divide 
the intervals [0, r) and [r, /3) into N and M sub-intervals, respectively, and use separate 
decompositions in each. We have 



e(r') 



E (T,,e*^'-^ , < r' < r 

r=-(N-l)/2 
{M-l)/2 

E Pre'"^^' , T <t' < (3 

I r=-(M-l)/2 



(104) 



with the reality condition a_.,. = a* , p_.r = p* and r-dependent frequencies Uj. = 2nr/T , z/^, 
27rr/(/3 — r). In terms of the new variables the representation ( |103| ) becomes 



V[a]V[p] exp 



X Tr ^ 



Texp 



X Texp 



- / dr' (k - xaoV - x E ^rC^^^^v) 
J dr' iK-xPoV-xT.Pre"''''^ 



(105) 



The discontinuous static-path approximation (DSPA), originally introduced in using 
the real-time framework, is now obtained by neglecting the contribution of the oscillations 
about the static paths (|102| ) which results in the two-dimensional integral 
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g(DSPA)^. 



1 



X 



ZiDSPA) 

xTr 



The partition function in this approximation in given by 



^(DSPA) 



x_ 

27r 



g-(/3-r)(A--XPoV)g- 



•T(/^-Xf70V) 



(106) 



(107) 



Note that Z'^^^^^^ acquires a r-dependence. 

In order to carry out the imaginary-time technique of Fourier-transforming Q (r) to get 
Q{iuJn) and extract the strength function G{uj) by an analytic continuation as was done 
above, it is necessary to obtain the functional dependence on r in (|106|) analytically. How- 
ever, the r-dependence of the traces involved is non-trivial and had to be studied numerically 
in Ref. even for the simple U{2) model; note that the analogous situation in the real-time 
framework (namely that the t-dependence of G^^^^^\t) is not given analytically) does not 
pose a problem since a numerical Fourier transform produces G^^^^^\uj) directly. Further- 
more, unlike the SPA case where we had a static Hamiltonian ho and Wick's theorem could 
be used to calculate the traces, here we have a discontinuous Hamiltonian 



ho{r') = 



K - xaoV , 
K - XPoV , 




(108) 



for which Wick's theorem is not applicable. In particular, without a generalization of Wick's 
theorem to this situation, it would be difficult to use our methods to calculate corrections 
due to small oscillations about the discontinuous paths. 



VI. CONCLUSION 

In this paper we present an approximation scheme, the PSPA, for the calculation of 
thermodjTiamic quantities and finite temperature response functions in finite fermionic sys- 
tems. The approximation is derived in the framework of the auxiliary-field path integral. 
We use an imaginary-time formulation which facilitates the extension of this approximation 
to physical quantities beyond the free energy and the level density to which it was previously 
limited. 

Testing the PSPA in a simple many-body model, we find that it improves on the SPA 
and is a good approximation for expectation values of observables as well as for low moments 
of strength functions. This indicates that the contribution of time-dependent fiuctuations 
about the static fields (neglected in the SPA) is significant. The required computational work 
involved in the PSPA includes a g-dimensional numerical integration over the static fields ctq 
and a diagonalization of a g x g matrix at each ctq -point, where q is the number of separable 
interactions in the Hamiltonian. This approximation breaks down at low temperatures when 
the small-oscillation correction factor diverges for the dominant static fields, indicating that 
large time- dependent fiuctuations become important. However, the breakdown occurs at 
temperatures well below the mean-field transition and does not affect the usefulness of the 
PSPA except at very low temperatures. For the strength function itself the PSPA results 
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become less reliable when the contribution of static paths with large discontinuity (at the 
time where the response function is calculated) is important. Further improvement would 
require the inclusion of both discontinuous static paths and small time-dependent oscillations 
around them. 

It would be interesting to test the PSPA methods for more realistic nuclear interactions, 
such as pairing plus multipole interactions. For these interactions, the SPA works better 
in a mixed pairing-density decomposition than in a pure density decomposition. Thus, it 
would be useful to extend the present PSPA techniques to such a mixed pairing-density 
decomposition in the HS representation. 

This work was supported in part by the Department of Energy Grant DE-FG02- 
91ER40608. 



APPENDIX: RPA FREQUENCIES AT FINITE TEMPERATURE 



The finite temperature RPA equations for a general two-body interaction Uijki are given 
by i 



ij } 



(Al) 



k,l 



where Aj^ = ej — ej and fi are the Fermi-Dirac occupation numbers /j = (1 + e^^')~^. The 
solutions of (|A1|) are the RPA frequencies Qi, and C,^ are the associated RPA amplitudes. 
The single-particle energies in ( |A1| ) correspond to the mean-field solution o"o but in the 
following we replace ctq by a general static field ctq. 

For a separable interaction as in with V Hermitean, or more generally for an interac- 
tion as in which is a sum of q such separable terms, Uujk = — J2 XaVfjVfj., and Eq. ( |A1|) 

can be rewritten as 



kl 



kl 



2]) can be converted to a set of coupled equations for the q quantities 

Ca^T.<kifk-mkl 



kl 



by multiplying Eq. (|A3|) by — fj) and summing over ij for each a' = 1, 

fi ~ fj 



E 



Sa'a + Xa E '"jl'^Tj 



An + 



C = 0; a'=l,...,q. 



(A2) 

(A3) 
, q. We get 
(A4) 



For ([A4| ) to have a non-trivial solution (i.e. not all Ca = 0), we require that the determinant 
of the coefficient matrix vanishes 



det 



r I a' a fi fj 



(A5) 
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Regarding the l.h.s. of Eq. (^) as a function of c<j, where we have substituted u for VL^, 
we notice that its roots are ±fi,y, while its poles are ±Ajj. It then follows that 



det 



r I a' a fi f 

Oa'a + Xa 2^ ^ji V. 



ml - u^') 
n '(A?. - 



(A6) 



Upon the substitution u; — > ia; we find 



det 



a fi fj 

Aij + iuj 



mi + 
n '( A?. + 



(A7) 



Using Eq. ( |A7|) for uj = ujr one obtains ( |36D or (|4l|) when one or several separable terms in 
the interaction are present, respectively. 



APPENDIX: THE MODEL 

The formalism developed in this paper is illustrated and tested in a simple Fermionic 
system, a variant of a model introduced in which is based on a U{2) algebra and is 



therefore solvable. This is a two-level system where each level is (7-fold degenerate and may 
therefore contain between zero and 2g Fermions. The Hamiltonian is given in terms of 
quasi-angular momentum operators 

H = 2eJ, - 2xJl (Bl) 

which has the form (|I]) with K = 2eJz and V = 2Jx- The quasi-angular momentum operators 
are given by 

9 , , ^9 



1 % 

i=l i=l 
\ 9 9 

= nJ2 - 4i«2i) , N = J2 {aliaii + alia2i) (B2) 



i=l i=l 



with N being the particle-number operator. The 2^^ states are arranged in f/(2)-multiplets 
{n,j) where the quantum numbers are the number of particles n = 0, 2g and the quasi- 
angular momentum j = 0, ...,n/2 or 1/2, ...,n/2 (depending on whether n is even or odd). 
Each multiplet (n, j) contains 2j + 1 states labeled by | njm) with m = are the 

eigenvalues of J^. To find the number of (n, j)-multiplets we first observe that the 

number of states with given n and m is 

^"(")=(n/2-»)(n/2 + ,J • f^^) 
Since every multiplet {n,j') with j' > j contributes a single state with = j, we also have 



27 



Therefore 



n/2 
j'=j 



(B4) 



n/2-jJ \n/2+j J \n/2 - j - ij \n/2 + j + 1 



(B5) 



which checks to give 



E (2i + lK(i) 

j=0(l/2) 



2^ 
n 



(B6) 



the total number of states with n particles. 

The Hamiltonian matrix in this basis is block-diagonal with identical blocks of 

dimension 2j + 1 for each pair {n,j), whose diagonalization gives the energies Ej„i and 
corresponding eigenstates. The exact partition function in the grand canonical ensemble is 
then 



2g n/2 j 
n=0j=0(l/2) m=-j 

The expectation value of an operator O is given by 

1 2g n/2 j 

{0) = ^i: E dnU) E {J^ I O I jm)e-'^(^--'^-) . 

^ n=0 j=0(l/2) m=-j 

The strength function associated with an operator O is given by 

1 2g n/2 j 

G{oj)^jYl E ^n(j) E \ ijm' \ O \ jm) \' e-^^'^^r.-.n) 

n=0 j=0(l/2) mm'=-j 

X ^ Im 



(B7) 



(B8) 



1 - e-l^^ 



(B9) 



which reduces to a sum over 5-functions as 77 0"*". 
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FIGURES 




FIG. 1. Free energy F{I3) = —f3 ^ log Z as a function of (3 for different values of k (see (^7[)) 
The SPA (dotted) and PSPA (dashed) results are obtained using (||) in Eqs. (|6|) and (|3|). Th( 
MFA result (dashed-dotted) is given by a steepest-descent treatment of the SPA integral. Thf 
exact result (solid) is calculated from (B7). 
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FIG. 2. Effective static-field free energy Fq{(3; cjo) (^) as a function of do in the SPA (dotted) 
and PSPA (dashed) at different temperatures /?. Shown is the case k = 1.5 where the mean-field 
phase-transition occurs at Pc = 1-61. 
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FIG. 3. Expectation value of as a function of /3 for different values of k. The SPA result 
(dotted) is given by ( pSj ) and the MFA result (dashed-dotted) is obtained from it by steepest descent. 
The PSPA result (dashed) is calculated using (|4^)-(|45|) and (|62|)- (|6^). These approximations are 
compared with the exact result ( |B8D (solid). 
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/c=1.5 



/c=3.0 




FIG. 4. Expectation value of J^, Jy and as functions of /? for different values of k. The 
SPA result (dotted) is given by (|7^ and the MFA result (dashed-dotted) is obtained from it by 
steepest descent. The PSPA result (dashed) is calculated using ([431)- (^5]) and (p8|)-(|7lD. These 
approximations are compared with the exact result ( |B8| ) (solid). 
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FIG. 5. Expectation value of the Hamiltonian H = 2eJz — '^xJx as a function of (3 for 

different values of k. Shown are the SPA (dotted), MFA (dashed-dotted), PSPA (dashed) and 
exact (solid) results. 
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oo 

FIG. 6. First moment of the strength function Mi = J dco u;G{uj) as a function of (3 for 

— oo 

different values of k, using ( |101D with D = J^- Shown are the SPA (dotted), MFA (dashed-dotted), 
PSPA (dashed) and exact (sohd) results. 
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/c= 1.5 



/c=3.0 




FIG. 7. Strength function G{uj) for D = for different values of k and at temperatures 
above (top), at (middle) and below (bottom) the mean- field transition (which does not occur for 
K = 0.5). Results are obtained using (|8|)-®) and (|2|)-(|6|) with t] = 0.1. Shown are the SPA 
(dotted), MFA (dashed-dotted), PSPA (dashed) and exact (^) (sohd) results. 
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FIG. 8. Same as in Fig. but for D = 
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